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The accurate but expensive product of geminals ansatz may be approximated by a geminal power, 
but this approach sacrifices size consistency. Here we show both analytically and numerically that 
a size consistent form very similar to the product of geminals can be recovered using a network 
of location specific Jastrow factors. Upon variational energy minimization, the network creates 
particle number projections that remove the charge fluctuations responsible for size inconsistency. 
This polynomial cost approach captures strong many-electron correlations, giving a maximum error 
of just 1.8 kcal/mol during the double-bond dissociation of H2O in an STO-3G atomic orbital basis. 



The overwhelming majority of electronic structure 
methods applied today rely fundamentally on the inde- 
pendent particle approximation (IPA). These methods, 
which include density functional theory [JJ , coupled clus- 
ter theory [2J, configuration interaction [3J, and many 
body perturbation theory [3J, all assume that the wave 
function is well approximated by a single Slater determi- 
nant (SD) in which the only correlations between elec- 
trons are those due to Fermi statistics. This assumption 
fails dramatically in a number of important cases dis- 
playing strong correlation between electrons, including 
multiple-bond breaking, excited states, transition metal 
compounds, and lattice Hamiltonians used in the study 
of high temperature superconductivity. While this failure 
can in some cases be rectified by active space methods 
that employ linear combinations of determinants, these 
methods' costs increase exponentially with system size. 
Indeed, when developing methods to treat strong corre- 
lation, one prefers to retain the formal properties of the 
SD: polynomially scaling cost, energies that are varia- 
tional (i.e. upper bounds), and size consistency, in which 
two non-interacting systems give the same total energy 
when modeled separately or together. 

One approach to this ideal is to generalize the SD, 
which is a product of one-particle functions (orbitals), to 
a product of two-particle functions (geminals) , known as 
the antisymmetric product of geminals (APG). 
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|*APG> = II &l°> & =X>«4ta!j. W 
i=l rs 

Here each operator Gi creates a pair of opposite-spin elec- 
trons in a two-particle geminal defined by the weights g % rs 
and operators and a/ s , that create f and 4- electrons 
in the sites (or orbitals) r and s. (The conclusions in this 
Letter generalize to same-spin pairs and pfaffians [U [5J , 
but to avoid unnecessary complication we restrict our- 
selves to opposite-spin pairs.) If no restrictions are placed 
on the form of the geminals, the resulting wave function 
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has been shown to be highly accurate. [SJ [7] However, 
the author is not aware of any polynomial cost, varia- 
tional methods for working with the general APG, and 
indeed it is more often approximated by requiring that 
the geminals be built from separate, mutually orthogo- 
nal sets of one-particle functions (APSG) [8HTU]. result- 
ing in methods such as perfect pairing (PP) [TT] and the 
resonating valence bond (RVB) [12]. While APSG meth- 
ods can achieve size consistency, variational energies, and 
polynomial cost, they lack correlation between electron 
pairs |13j and are thus unsuitable for treating strong cor- 
relations between more than two electrons [2]. While 
corrections can be applied via configuration interaction 
[TOl fT5] , coupled cluster [IH [TBI Hi] > perturbation theory 
[TU] , and Hopf algebra [TS] [T5] , none of these approaches 
simultaneously retain polynomial cost, variational ener- 
gies, and size consistency. 

Building on the work of Casula and Sorella (see Refs. 
[2OH22] and especially [23]), we present an ansatz that 
captures strong inter-pair correlations while retaining 
polynomial cost, variational energies, and size consis- 
tency. To the best of our knowledge, this is the first 
example of a method that achieves all of these properties 
for a general system and an ab initio Hamiltonian. 

Ansatz. — We begin our construction with the well- 
known geminal power (AGP) ansatz, 

|*agp> = ^ /2 |0) F = J2 frsal t al ~ ]T G h (2) 

rs i 

in which the (bosonic) electron pairs all reside in the same 
low-energy geminal F, which should be similar to the sum 
of the APG geminals Gi. For those more familiar with 
the superconducting Bardeen-Cooper-Schrieffer (BCS) 
ansatz [24 , it may be helpful to consider that the AGP is 
the RVB equivalent of a particlc-number-projected BCS, 
with the real space pairing matrix f rs related by a Fourier 
transform to the BCS fc-space weights (see Ref. [25], Eqs. 
4.9-4.10). While the AGP admits a number of polyno- 
mial cost, variational methods [2U1 HU I25H3TJ] (one of 
which 30J achieves a mean-field n 3 cost), it suffers from 
a severe size consistency problem resulting from terms in 
which a single operator Gi is repeated, placing four or 
more electrons in the same local geminal. 
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FIG. 1: (color online) A cartoon schematic of our JAGP 
ansatz, in which the geminal power is constructed from non- 
orthogonal local geminals (ovals) describing bonds between 
neighboring atoms (circles). Particle number projectors (rect- 
angles) built from Jastrow factors constrain electron counts 
on atoms and groups of atoms to remove the "ionic" AGP 
terms responsible for size consistency errors. 



Our approach is to eliminate these charge transfer or 
"ionic" terms by enforcing local particle number distribu- 
tions with a network of Jastrow factors. We thus produce 
a size consistent Jastrow- AGP (JAGP) ansatz, 



|* JAGP > = exp ( E -W ] ^ /2 |°>' ( 3 ) 

pq (Tretd 

in which the Jastrows operate on the bare AGP in the 
same way they operate on a bare SD in the traditional 
Jastrow-Slater ansatz [31] . The Jastrow factors J Va q T 
inspect the occupation (00, 01, 10, or 11) of each orbital 
pair p a q r and apply a corresponding scalar factor to the 
wave function to favor or penalize the configuration of 
that particular pair. They are defined by 
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where C is a tensor of penalty factors and the operator 
PP" gives one if orbital p a has occupation n and zero 
otherwise. Note that together, these Jastrow factors are 
equivalent to the correlator product state tensor network 
[3"2l 155] , and so we refer to them as a location-specific 
Jastrow factor network. In this Letter, we demonstrate 
that this network restores size consistency to the AGP 
by imposing local particle number constraints, delivering 
an ansatz that is similar in character to the APG, size 
consistent, variational, polynomial cost, and effective at 
treating strong many-electron correlations. 

Charge fluctuations. — To be size consistent, a wave 
function must factor into a product of subsystem wave 
functions \^ab) = l^yOl^s) when applied to two non- 
interacting subsystems A and B. As noted previously by 
Sorella, Casula, and Rocca [53] , unphysical charge fluctu- 
ations prevent this factorization. For a simple example, 
imagine two H2 molecules described by geminals G a and 
G B - The AGP built from these geminals, (G A + G B ) 2 \0), 
contains both the neutral term GaGb |0) and the unphys- 
ical ionic terms G^JO) and G B \0) in which all four elec- 
trons reside on one molecule. Without the ionic terms, 



this AGP would factor correctly and be size consistent. 

We may generalize this analysis by expanding the AGP 
in the basis of occupation number vectors |n), each of 
which specifies a unique occupation pattern of the or- 
bitals. 
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Here the coefficients simplify to determinants of the oc- 
cupied pairing matrices ^ n [34], which are obtained by 
deleting from / rows and columns corresponding to un- 
occupied orbitals. Note that the sum is restricted to 
states with the correct total t an d I electron counts 
N t = N l = N/2. 

An intuitive guess for \^ab) is to take the AGP gem- 
inal as the sum of the subsystem geminals and the Jas- 
trow factor as the product of the subsystem Jastrows, in 
which case the pairing matrix / will be block diagonal 
with blocks equal to the subsystem matrices /a and fs , 
and the Jastrows will be defined by C = Ca + Cb- Such 
a choice results in 
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where the determinant factors due to the block- 
diagonality of / and the Jastrow factors due to the ad- 
ditive separability of C. However, I^ab) does not factor 
into I VPyi) I SE'^b) , because the summation over orbital oc- 
cupations n contains ionic terms in which electrons are 
transferred between subsystems. 

Using real space three-body Jastrow factors, Sorella et 
al showed [23] that these spurious charge fluctuations can 
be partially suppressed, mitigating the size consistency 
error. However, removing the error completely through 
this approach would require unlimited flexibility in the 
Jastrow. In practice, their wave function retained a size 
consistency error on the order of leV in the carbon dimer 
[23], although the effect on binding energies was much 
smaller due to error cancellation. We expand on this 
idea, showing how Jastrow factors can eliminate the size 
consistency error entirely. 
Partial number projection. — Consider the operator 



Q(a, M, X) = exp 
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which we call a partial number projection operator fa- 
voring M electrons in the set of orbitals X. In the limit 
a — > 00, this becomes a strict projection, deleting terms 
in which X's electron count differs from M. We may thus 
fix the subsystem electron counts and delete ionic terms 
using the operators Qa — Q(ce, Na^, A^)Q(a, Na±, A±) 
and Qb = Q(a, Nb^-, B^)Q(a, Nb^, B^), which when ap- 
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plied to Eq. ([6| produce the desired factorization. 

lim Q a Qb\*ab) = \*a)\*b) (8) 
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Thus if we can apply appropriate projections, JAGP will 
factor and be size consistent. 

The advantage of our ansatz is that the partial projec- 
tion operators Qa and Qb can be built into the Jastrow 
network. To see how, expand the square in Eq. ^ and 
drop the constant term exp(— aM 2 ), which only affects 
normalization, to obtain 



Q(a,M,X) 
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exp 
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where (3 = 2Ma/k and k is the number of orbitals in 
X. Inspecting Eq. ( |To| reveals that the Jastrow network 
defined in Eqs. ([3| and Q can contain any combination 
of partial projection operators. The JAGP is therefore 
capable of deleting ionic terms by restricting subsystem 
electron counts, making it factorizable and size consis- 
tent. Furthermore, if we take our AGP geminal as a sum 
of the localized but non-orthogonal APG geminals, par- 
tial number projections can help ensure that each local 
geminal has the correct number of electrons. Our JAGP 
thereby emulates the structure of the APG. 

Variational minimization. — We use variational Monte 
Carlo (VMC) [3"Tll3"5] to evaluate and minimize the JAGP 
energy by varying independently all elements of the pair- 
ing matrix / and Jastrow factor penalty tensor C. The 
Hamiltonian is the typical Born-Oppcnheimer approx- 
imation to the electronic Hamiltonian with relativistic 
terms neglected. Note especially that we work in Fock 
space rather than real space. We use an improved ver- 
sion of the Linear Method optimizer along the lines pro- 
posed in Ref. [36], the details of which will be presented 
elsewhere [37] . For the present discussion, it suffices to 
convey that this method is variational with a cost of 
Oinsn^n^), where n s , n Q , and n u are the sample size 
and the numbers of occupied and unoccupied orbitals. 

Hydrogen gas. — A collection of n well separated hydro- 
gen molecules reveals the severity of AGP's charge fluctu- 
ations. Working in a symmetrically orthogonalized STO- 
3G basis [3H], in which a single Is orbital is centered on 
each H, we may define the AGP geminal as a sum of PP 
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FIG. 2: (color online) The average number of unphysical 
charge transfers per molecule in a system of n well separated 
H2 molecules. The wave function is a PP-parameterized AGP 
with various partial number projections. The dotted line is a 
fit showing the asymptotic 1/n decay for a = 0.1. 
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FIG. 3: (color online) Energy errors per molecule for n well- 
separated H2 molecules. For AGP, both the PP and optimized 
versions of the wave function are shown. 



geminals, 
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Here x 2 + y 2 = 1, Q is a partial number projection oper- 



ator suppressing charge fluctuations, and g\^n an( l u i\/i 
create electrons in the bonding and antibonding orbitals, 
respectively, of the ith H2 molecule. If we parameterize 
Q to apply a penalty of e~ 2a for each incorrect H 2 elec- 
tron count, then the average number of charge transfers 
(i.e. the number of [H 2 ] 2+ ions) will be 
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FIG. 4: Energy errors relative to FCI for the symmetric 
dissociation of minimal basis H2O with bond angle 109.57°. 
JAGP's statistical uncertainties are smaller than the symbols. 



where the contributions are grouped by the number of 
charge transfers I. Figure [2] shows that (Ncr/n) decays 
as X/n in the thermodynamic limit n — > 00, recovering 
the size extensivity of the BCS ansatz. However, the 
steep growth of (Ncr/n) for small n is unacceptable for 
quantum chemistry, where system sizes range from tens 
to hundreds of bonding electron pairs. Encouragingly, 
(Ncr/n) is very sensitive to increasing a, showing that 
charge fluctuations are easily suppressed in our ansatz. 

Without partial number projection, the charge fluctu- 
ations in an AGP built from PP geminals render it less 
accurate than the IPA for small n, as shown in Figure 
[3] Here we use the somewhat more realistic symmet- 
rically orthogonalized 6-3 1G basis [33] • The errors for 
the variationally optimized AGP are less than those of 
the IPA, but they remain large compared to those of sin- 
gles and doubles configuration interaction (CISD), whose 
well-known size consistency problem is much less severe. 
Most importantly, variational optimization (with initial 
guess / = random, C = 0) of our JAGP completely 
removes size consistency errors and produces the exact 
PP result. This shows our optimization can discover the 
need for particle number projection and impose it auto- 
matically. 

It is worth noting that a significant component of the 
JAGP's correlation energy is size extensive (i.e. it scales 
linearly with system size for large systems), because the 
JAGP can always be made to contain PP, and PP ener- 
gies are size extensive. Less clear is whether the entire 
JAGP energy is extensive, which clearly merits further 
investigation. 

Double bond dissociation. — To demonstrate JAGP's 
ability to capture strong inter-pair correlations while 
maintaining size consistency, we have applied it to the 
symmetric double-bond dissociation of H 2 in a symmet- 
rically orthogonalized STO-3G basis. We first optimized 
the wave function for a single molecule, starting from a 



TABLE I: JAGP energies for collections of n well separated 
H2O molecules with bond lengths 1.4A and angles 109.57°. 
Statistical uncertainty in final digit given in parentheses. 



n 


E/n (a.u.) 


1 


-74.90371(1) 


2 


-74.90374(3) 


4 


-74.90369(3) 


8 


-74.90376(5) 



very poor initial guess (/ = random, C = 0). Figure [4] 
shows that the maximum error relative to full configu- 
ration interaction (FCI) is 1.8 kcal/mol, a factor of 2.5 
smaller than the 4.5 kcal/mol error produced by unre- 
stricted coupled cluster with singles, doubles, and per- 
turbative triples (UCCSD(T)). In terms of correlation 
energies (defined with respect to an unrestricted SD), 
JAGP retains above 90% of the correlation across the 
whole curve, while UCCSD(T) is less balanced with cor- 
relation recovery ranging from over 99% near equilibrium 
down to 75% upon dissociation. 

After optimizing our ansatz for one water molecule, we 
tested size consistency by constructing wave functions for 
two, four, and eight well separated water molecules. The 
overall geminals were built as sums of monomer geminals, 
and the Jastrow tensor C as the sum of the monomers' 
plus the terms necessary to impose partial number pro- 
jection with a = 2 on the f an d 4 electron occupa- 
tions of each molecule. Table |T] reveals that the energy 
per molecule is independent of the number of molecules, 
showing that JAGP is size consistent even when it is not 
exact (as was the case for H2). Finally, note that the 
eight water case corresponds to a 40-orbital active space. 

Conclusions. — We have shown that a geminal power 
augmented with a network of location-specific Jastrow 
factors recovers size consistency in a localized one parti- 
cle basis, producing an ansatz similar in character to the 
powerful but expensive product of non-orthogonal gem- 
inals. The resulting method is variational, size consis- 
tent, polynomial cost, and able to capture strong many- 
electron correlations. It completely removes unphysical 
charge fluctuations from a dilute H2 gas and accurately 
captures the strong correlations of water's double-bond 
dissociation. We believe it is the first geminal method 
satisfying all of these properties and that it is a promis- 
ing candidate for applications to other strongly correlated 
systems. 
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